function y = func_delta(p, m, CONSTS)
    
    k0    = CONSTS.k0;
    a     = CONSTS.a;
    eta   = CONSTS.eta;
    eps_a = CONSTS.eps_a;

    s      = func_s(p, CONSTS);
    n1     = func_n1_from_p(p, CONSTS);
    n2     = func_n2_from_p(p, CONSTS);
    K_m    = func_K_m(p, m, CONSTS);
    J_m1   = func_J_m1(p, m, CONSTS);
    J_m2   = func_J_m2(p, m, CONSTS);
    J_m1_1 = func_J_m1_1(p, m, CONSTS);
    J_m2_1 = func_J_m2_1(p, m, CONSTS);

    if (m>0)
        y = (n2-n1).*((p.^2)./eps_a).*((m^2)./((k0*a*s).^4)) - (n2-n1).*(K_m).^2 + ...
            n1.*(eta/eps_a).*J_m1_1.*J_m2_1 - n2.*(eta/eps_a).*J_m1.*J_m2 - ...
            n2.*J_m1.*K_m + (eta/eps_a)*n1.*J_m1_1.*K_m + n1.*J_m2_1.*K_m - ...
            (eta/eps_a)*n2.*J_m2.*K_m + p.*(eta/eps_a).*(m./((k0*a*s).^2)).* ...
            (-J_m2_1 + J_m1 + J_m1_1 - J_m2);
    elseif(m<0)
        m1=abs(m);
        y = (n2-n1).*((p.^2)./eps_a).*((m1^2)./((k0*a*s).^4)) - (n2-n1).*(K_m).^2 + ...
            n1.*(eta/eps_a).*J_m1_1.*J_m2_1 - n2.*(eta/eps_a).*J_m1.*J_m2 - ...
            n2.*J_m1.*K_m + (eta/eps_a)*n1.*J_m1_1.*K_m + n1.*J_m2_1.*K_m - ...
            (eta/eps_a)*n2.*J_m2.*K_m - p.*(eta/eps_a).*(m1./((k0*a*s).^2)).* ...
            (-J_m2_1 + J_m1 + J_m1_1 - J_m2);
    else
        y = (n2-n1).*((p.^2)./eps_a).*((m^2)./((k0*a*s).^4)) - (n2-n1).*(K_m).^2 + ...
            n1.*(eta/eps_a).*J_m1_1.*J_m2_1 - n2.*(eta/eps_a).*J_m1.*J_m2 - ...
            n2.*J_m1.*K_m + (eta/eps_a)*n1.*J_m1_1.*K_m + n1.*J_m2_1.*K_m - ...
            (eta/eps_a)*n2.*J_m2.*K_m + p.*(eta/eps_a).*(m./((k0*a*s).^2)).* ...
            (-J_m2_1 + J_m1 + J_m1_1 - J_m2);
    end
        

end

function y = func_J_m1(p, m, CONSTS)
    k0 = CONSTS.k0;
    a = CONSTS.a;
    q1 = func_q1_from_p(p, CONSTS);
    alfa_1 = func_alfa_1(p, CONSTS);

    if(m>=0)
        bessel_m_q1 = besselmx(double('J'), m, k0*a*q1, 0);
        bessel_mp1_q1 = besselmx(double('J'), m+1, k0*a*q1, 0);

        y = (bessel_mp1_q1./((k0*a*q1).*bessel_m_q1)) + ...
               ((alfa_1.*m)./((k0*a*q1).^2));
    else
        m1=abs(m);
        bessel_m_q1 = besselmx(double('J'), m1, k0*a*q1, 0);
        bessel_mp1_q1 = besselmx(double('J'), m1+1, k0*a*q1, 0);

        y = (bessel_mp1_q1./((k0*a*q1).*bessel_m_q1)) - ...
               m1*((alfa_1+2)./((k0*a*q1).^2));
    end

end

function y = func_J_m2(p, m, CONSTS)
    k0 = CONSTS.k0;
    a = CONSTS.a;
    q2 = func_q2_from_p(p, CONSTS);
    beta_2 = func_beta_2(p, CONSTS);

    if(m>=0)
        bessel_m_q2 = besselmx(double('J'), m, k0*a*q2, 0);
        bessel_mp1_q2 = besselmx(double('J'), m+1, k0*a*q2, 0);

        y = (bessel_mp1_q2./((k0*a*q2).*bessel_m_q2)) - ...
               ((beta_2.*m)./((k0*a*q2).^2));
    else
        m1=abs(m);
        bessel_m_q2 = besselmx(double('J'), m1, k0*a*q2, 0);
        bessel_mp1_q2 = besselmx(double('J'), m1+1, k0*a*q2, 0);

        y = (bessel_mp1_q2./((k0*a*q2).*bessel_m_q2)) + ...
               m1*((beta_2-2)./((k0*a*q2).^2));
    end
   
end

function y = func_J_m1_1(p, m, CONSTS)
    k0 = CONSTS.k0;
    a = CONSTS.a;
    q1 = func_q1_from_p(p, CONSTS);
    beta_1 = func_beta_1(p, CONSTS);

    if(m>=0)
        bessel_m_q1 = besselmx(double('J'), m, k0*a*q1, 0);
        bessel_mp1_q1 = besselmx(double('J'), m+1, k0*a*q1, 0);

        y = (bessel_mp1_q1./((k0*a*q1).*bessel_m_q1)) - ...
               ((beta_1.*m)./((k0*a*q1).^2));
    else
        m1=abs(m);
        bessel_m_q1 = besselmx(double('J'), m1, k0*a*q1, 0);
        bessel_mp1_q1 = besselmx(double('J'), m1+1, k0*a*q1, 0);

        y = (bessel_mp1_q1./((k0*a*q1).*bessel_m_q1)) + ...
               m1*((beta_1-2)./((k0*a*q1).^2));
    end
end

function y = func_J_m2_1(p, m, CONSTS)
    k0 = CONSTS.k0;
    a = CONSTS.a;
    q2 = func_q2_from_p(p, CONSTS);
    alfa_2 = func_alfa_2(p, CONSTS);

    if(m>=0)
        bessel_m_q2 = besselmx(double('J'), m, k0*a*q2, 0);
        bessel_mp1_q2 = besselmx(double('J'), m+1, k0*a*q2, 0);

        y = (bessel_mp1_q2./((k0*a*q2).*bessel_m_q2)) + ...
               ((alfa_2.*m)./((k0*a*q2).^2));
    else
        m1=abs(m);
        bessel_m_q2 = besselmx(double('J'), m1, k0*a*q2, 0);
        bessel_mp1_q2 = besselmx(double('J'), m1+1, k0*a*q2, 0);

        y = (bessel_mp1_q2./((k0*a*q2).*bessel_m_q2)) - ...
               m1*((alfa_2+2)./((k0*a*q2).^2));
    end
end

function y = func_K_m(p, m, CONSTS)
    k0 = CONSTS.k0;
    a = CONSTS.a;
    s = func_s(p, CONSTS);

    if(m>0)
        y = (besselk(m+1, k0*a*s, 0)./((k0*a*s).*besselk(m, k0*a*s, 0))) - ...
               (m./((k0*a*s).^2));
    elseif(m<0)
        m1=abs(m);
        y = (besselk(m1+1, k0*a*s, 0)./((k0*a*s).*besselk(m1, k0*a*s, 0))) - ...
               (m1./((k0*a*s).^2));
    else
        y = (besselk(m+1, k0*a*s, 0)./((k0*a*s).*besselk(m, k0*a*s, 0)));
    end
end
